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Abstract 

This paper analyzes the random fluctuations obtained by a heterogeneous multi- 
scale first-order finite clement method applied to solve elliptic equations with a random 
potential. Several multi-scale numerical algorithms have been shown to correctly cap- 
ture the homogenized limit of solutions of elliptic equations with coefficients modeled 
as stationary and ergodic random fields. Because theoretical results are available in 
the continuum setting for such equations, we consider here the case of a second-order 
elliptic equations with random potential in two dimensions of space. 

We show that the random fluctuations of such solutions are correctly estimated by 
the heterogeneous multi-scale algorithm when appropriate fine-scale problems are solved 
on subsets that cover the whole computational domain. However, when the fine-scale 
problems are solved over patches that do not cover the entire domain, the random 
fluctuations may or may not be estimated accurately. In the case of random potentials 
with short-range interactions, the variance of the random fluctuations is amplified as 
the inverse of the fraction of the medium covered by the patches. In the case of random 
potentials with long-range interactions, however, such an amplification docs not occur 
and random fluctuations are correctly captured independent of the (macroscopic) size 
of the patches. 

These results are consistent with those obtained in [8] for more general equations in 
the one-dimensional setting and provide indications on the loss in accuracy that results 
from using coarser, and hence computationally less intensive, algorithms. 

Keywords: Equations with random coefficients, multi-scale finite element method, 
heterogeneous multi-scale method, corrector test, long-range correlations. 

AMS subject classification (2010): 35R60, 65N30, 65C99 

1 Introduction 

Differential equations with highly oscillatory coefficients arise naturally in many areas of 
applied sciences. The microscopic details of such equations are difficult to compute. Nev- 
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ertheless, when the heterogeneous medium has certain properties involving separation of 
scales, periodicity, or stationary ergodicity, homogenization theories have been developed 
and they provide macroscopic models for the heterogeneous equations; see e.g. [19 ] l2i j [25], 
Many multi-scale algorithms have been devised to capture as much of the microscopic scale 
as possible without solving all the details of the micro-structure [21 [T5l [Til I18j . Such a 
scheme is viewed as correct if it can well approximate the macroscopic solution when the 
heterogeneous medium satisfies conditions for homogenization to happen. Homogenization 
theory thus serves as a benchmark which ensures that the multi-scale scheme performs well 
in controlled environments, with the hope that it will still perform well in non-controlled 
environments, for instance when ergodicity and stationarity assumptions are not valid. 

In many applications such as parameter estimation and uncertainty quantification, esti- 
mating the random fluctuations (finding the random corrector) in the solution is as impor- 
tant as finding its homogenized limit j9j[23]. When this is relevant, another benchmark for 
multi-scale numerical schemes that addresses the limiting stochasticity of the solutions is 
plausible: One computes the limiting (probability) distribution of the random fluctuation 
given by the multi-scale algorithm in the limit that the correlation length of the medium 
tends to while the discretization size h of the scheme is fixed. If this /i-dependent distri- 
bution converges, as h — > 0, to the limiting distribution of the corrector of the continuous 
equation (before discretization), we deduce that the multi-scale algorithm asymptotically 
correctly captures the randomness in the solution and passes the random corrector test. 

Such proposal requires a controlled environment in which the theory of correctors is 
available. We introduced and analyzed such a benchmark in [8] using an ODE model whose 
corrector theory was studied in [1 H [7] . The main purpose of this paper is to provide and 
analyze another benchmark using a PDE model whose corrector theory was studied in 
[5| I16[ IB]. This is necessary because many multi-scale schemes that are different in higher 
dimensions turn out to be equivalent for the aforementioned ODE model, e.g. those in [TJ 
and [18]. In the rest of this introduction, we first review some main results in [8]. Then we 
introduce the results of the current paper that address the corrector test using an elliptic 
PDE with random potential. 

1.1 Corrector test using an ODE with random elliptic coefficient 

The corrector test is based on the homogenization and corrector theory of the following 
equation: 



Here, the diffusion coefficient a(j,iv) is obtained by rescaling a(x,co) which is a random 
process on some probability space (O, P). It is well known |21| I25] that (and this general- 
izes to higher dimensions as well) when a(x,co) is stationary, ergodic, and uniformly elliptic, 
then the solution u e converges to the following homogenized equation with deterministic 
and constant coefficient: 
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In the one-dimensional case, the coefficient a* is the harmonic mean of a(x), i.e., the inverse 
of the expectation of a -1 . We denote by q(x) the deviation of l/a(x) from its mean 1/a*. 
The corrector theories for the limiting distribution of u £ — u$ were studied by [TT] . 
The results in these papers are represented in path (Hi) of the diagram in Fig. [TJ The 
limiting distribution showing at the lower-right corner depends on the de-correlation rate of 
q(x). When q is strongly mixing with integrable mixing coefficient (see (|2.3p below), then 
13 = 1 and is a standard Brownian motion multiplied by a, a factor determined by the 
correlation function of q as detailed in (|2.2p below. When q has a heavy tail (is long-range) 
in the sense of (L1-L3) in section [21 we should take f3 = a, a < 1 being defined in (|2.4p . 
and W a is the fractional Brownian motion with Hurst index 1 — § multiplied by certain 
factor. These convergence results are understood as convergence in distribution in the space 
of continuous paths C([0, 1]). 

The corrector test for multi-scale numerical schemes is therefore the following: Let h be 
the discretization size and u £ (x) the solution to (jl.ip yielded by the scheme. Let Uq(x) be 
the solution yielded by the same scheme applied to (|1.2j) . The discrete corrector is u £ — Uq. 
According to the de-correlation property of q(x), we choose e 13 and interpret W 13 as before. 
We say that a numerical procedure is consistent with the corrector theory and that it passes 
the corrector test when the diagram in Fig. [T] commutes: 

Figure 1: A diagram describing the corrector test with a random ODE. 

«e - «q / \ h->0 U £ -Uq 
r^-{X,UJ) > 



s(> (i) Ve? 



e-s-0 



(ii) (iii) 



£->-0 



J L h (x,y)dW^(y) j \a*) 2 ^(x,y)^(y)dW^(y). 



More precisely, we need to characterize the intermediate limit in path (ii) which appears 
on the left of the diagram. In this step, h is fixed while the correlation length e is sent to 
zero. The intermediate limit distribution is /i-dependent. Very often, it can be described as 
a stochastic integral as shown and we need to determine the kernel function L h (x, y). Next, 
we need to verify the converge path (iv) which is taken as h — > 0. The numerical scheme is 
said to pass (or fail) the corrector test if this limit holds (or does not). 

In [8], we considered a Finite Element Method (FEM) based scheme in the framework 
of Heterogeneous Multiscale Methods (HMM), which is a general methodology for designing 
sublinear algorithms for multi-scale problems by exploiting special features of the problem, 
e.g. scale separation [15] . The macro-solver of this FEM-HMM scheme uses the standard 
PI element on a uniform grid of size h. The corresponding discrete bilinear form which 
approximates the continuous bilinear form associated to (jl.ip is 

A h (u h ,v h ) = Y^-(xi)a*^(x J )h^ f 1 ^-( x ) a *^(x)dx =: A(u\v h ). (1.3) 
^ dx J dx J Jn dx dx 

Here, a simple middle-point quadrature is used for the integral and Xj, j = 1, • • • , N = 1/h 
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are the evaluation points. Since the effective coefficient a* is unknown apriori, the FEM- 
HMM scheme approximates the discrete integrand by 

du h . , dv h . If du h dv h . . 

cte v ^ dx K v 5 Jr. x dx K ' y ' dx y ' 

where Ij$ = (xj — 5/2, Xj + 6/2) is a patch inside the discretization interval Ij = (xj — 
h/2, xj + h/2); the functions u h and v h are given in terms of {4^} where {(/>*} are the nodal 
bases and {(ft 1 } are given by the micro-solver 

;^ (x) ^ (x) = ' xeI ^ (i.4) 

(ftj(x) = <ftj{x), x £ dlj$. 

When 5 = h, this scheme coincides with those in |18|. Q]. It is known that one can choose 
5 < h to greatly reduce computational cost while still approximating the macroscopic 
solution quite well [15]. 

The main result of [8] shows that the corrector test for the above FEM-HMM scheme 
depends on the correlation structure of the random media. More precisely, for a long range 
correlated media (L1-L3 in section [2TT]) . the scheme is robust for the corrector test: the final 
limit in path (iv) of the diagram in FigQ] agrees with the theoretical Gaussian limit for all 
5 < h. For a short range correlated media (S1-S3 in section [2TTj) . however, this holds true 
only for 5 = h. The final limit for 5 < h is an amplified version of the theoretical Gaussian 
limit with an amplification factor (h/5) 1 ^ 2 , which shows that reducing the computational 
cost results in an amplification of the variance of the numerical calculations. 

1.2 Corrector test using elliptic PDE with random potential 

The main objective of this paper is to provide a two dimensional corrector test. Such a 
strategy generalizes to arbitrary space dimensions, although for concreteness, we concentrate 
on the two-dimensional setting. A full theory of random fluctuations for second order 
elliptic PDE with highly oscillating random diffusion coefficients in dimension higher than 
one remains open and we can not use it for the corrector test. Instead, we base the test on 
the following elliptic equation with random potential: 

f - Au £ + (q + q £ )u £ (x,uj) = f, x G Y, 
1 u E (x,u) = 0, x G dY. 

The coefficient in the potential term consists of a smooth varying function qo , and a highly 
oscillatory random function q(e~ l x,oj) denoted by q e (x) for simplicity. The random field 
q(x, oj) is assumed to be stationary ergodic and mean-zero. When e goes to zero, the solution 
u £ converges in L 2 (Q x Y) to the homogenized solution uq that solves 

f - Au + q u (x) = f, x G Y, 

\ u (x) = 0, x G dY. ^ ' ' 

The corrector theory for the above homogenization is well understood; see [I6l[5j[6]. When 
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Figure 2: A diagram describing the corrector test with a random PDE. 



U £ ~ Uq \ h^O 
U (x,U)),ip) ► 



' £ 



(<■<■) 



L h [4>]{ Xl y)dWP(y) 



h->-0 



(Hi) 



■(x,u)),(p 



(iv) Jy 



^(x)G(x,y)u (y)dW^(y). 



the corrector u £ — uq is properly scaled, it converges to a stochastic integral in a weak sense. 
This is described by the path (Hi) of the diagram in Fig. [2j Both the scaling factor and the 
limit depend on the correlation structure of the random field. These results are reviewed in 
Section [2] below. As in the ODE (one-dimensional) setting, a corrector test can be sketched 
as in the diagram of Fig. [2j For a given multi-scale scheme, which yields v!^ and Uq when 
it is applied to (|1.5p and (jl.6p . respectively, the main tasks are again to characterize the 
intermediate convergence in path (it) while e is sent to zero first and to check the validity 
of path (iv) when h is sent to zero afterwards. 



Figure 3: Left: Triangulation of the unit square. Right: Shrinking from K to Kg with 
respect to the barycenter. 




X\ X{ XJV x 



Now we introduce a heterogeneous multi-scale scheme for (|1.5p . The weak formulation 
of the equation is to find u £ in the Sobolev space Hq(Y) so that A £ (u £ ,v) = (f,v) for all 
v € Hq(Y). Here and below, (•, •) denotes the usual pairing; A £ is the bilinear form 

A £ (u,v) = J Vu-Vv + (q + q £ )uv dx, Vu, v £ Hq(Y). (1.7) 

Since we always assume that qo + q £ is positive, the weak formulation is well-posed thanks 
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to the Lax-Milgram lemma. The scheme that will be considered is based on FEM. For 
simplicity, Y is taken as the two dimensional unit square (0, l) 2 . Let Th be the standard 
uniform triangulation as illustrated in Fig. [3J Here, the typical length of the triangles is 
h = 1/N and N is the number of partitions on the axes. We consider first-order Lagrange 
elements. Associated to each (interior) nodal point (ih,jh), there is a continuous function 
(j^o which is linear polynomial restricted to each triangle K £ Th and which has value one 
at this nodal point and has value zero at all other nodal points. Note that the index i,j 
runs from 1 to N — 1. The space V h spanned by {</> Jjf } is a finite dimensional subspace of 
H^(Y). The hetero geneous multi-scale scheme for (jl.5p is to find u e ' € V h that satisfies 

A h /{u h /,v h ) = {f,v h ), for all v h £ V h , (1.8) 

where Ae' S is a bilinear form on V h x V h which approximates A e as follows: 

A h /{u h ,v h ) := Y \K\ (-^— I Vu h ■ Vv h + (q + q e )u h v h dx] . (1.9) 

Here, K$ C K is a patch centered at the barycenter of K and has typical length 5; the 
symbol | • | means taking the area. Ae' 5 hence is a numerical quadrature for the integral 
in (|1.7[) using averaged value around the barycenters of the elements. The scheme (jl.8p is 
analyzed in Section [3] and it is well-posed. 

When the above scheme is applied to the homogenized equation (jl.6p , it yields a solution 
Ug ,<5 in V h so that 



A^ 5 (4' 5 ,v h ) = (4' 5 ,v h ), for all v h 6 V h , (1.10) 



and Aq' S is given by 



A%\u h ,v h ) := y \K\ (JL [ Vu h ■ Vv h + q u h v h dx) . 

The discrete corrector function is defined to be the difference between u^ ,S and u^' S . 

Remark 1.1. The ratio is the two dimensional analog of 5/h in the aforementioned 

FEM-HMM scheme for the ODE setting. It measures savings in the computational cost. As 
in the ODE setting, we expect the corrector test to depend on the ratios {|i^j|/|.K"| : K £ 
Th}- To simplify notations, we assume ji^l/l-fTl is a constant for all K. For instance, 
consider a typical triangle K with vertices (0,0), (h,0) and (0, h). Kg can be obtained by 
shrinking with respect to the barycenter (h/3,h/3) so that it has vertices {{h — 5)/3, (h — 
5)/3),((h + 25) /3, (h - 5)/3) and ({h-5)/3,{h + 28) /3) ; see Fig. H If all K s are obtained in 
this manner, we have y^i^l/l-^l = {5/h) d / 2 with d = 2. More general patches than those of 
the paper could also be considered without changing our main conclusions; we concentrate 
on a specific construction of the small patches Kg to simplify notation and denote the ratio 
V#W by 5/h. 

1.3 Main Results 

The main results of this paper concern the limiting distribution of the discrete corrector 
Us ,S — Uq' 6 with proper scaling. They depend on the correlation structure of the random 
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field q £ . We refer to section [2~T1 below for notation. In particular, SRC (respectively LRC) 
stands for short (respectively long) range correlation. 

Theorem 1.2. Let u^ ,S and u^ ,S be the solutions obtained from the heterogeneous multi- 
scale schemes (jl.8p and (|1.10p . respectively. Assume that qo is a nonnegative constant and 
f is inC 2 (Y), i.e., twice continuously differentiable overY. For an arbitrary test function 
<p € C 2 (Y), the following holds. 

(1) In the SRC setting, as e goes to zero while h and 5 satisfying h > 5 S> e are fixed, we 
have 
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<p(x)[u*> s - u h /]dx dlst ™° n > a J L h [<p](x)dW(x). (1.11) 

Here L h [ip] is a bounded function defined in (|4,19p below. 

(2) As h and 5 go to zero with the ratio 5/h < 1 being fixed, we have 

a f L h [ip](x)dW(x) distributi ° n ) %a [ gcp(x)u (x)dW(x). (1.12) 
Jy ° Jy 

(3) In the LRC setting, the above convergences are modified by 

-±= [ v (x)[v^ ~ u h /]dx / L h M{x)W«{dx), (1.13) 

ye Jy e ^ u Jy 

and 

f L h [<p](x)W a (dx) distribution > v 7 ^ [ G<p{x)u (x)W a (dx). (1.14) 
Jy Jy 

Above, Q is the solution operator of (|1.6P ; n is defined after (12. 5p below. W is the standard 
multi-parameter Wiener process; W a (dy) is formally defined to be W a (y)dy and W a (y) is 
a Gaussian random field with covariance function given by ¥,{W a (x)W a (y)} = \x — y\~ a . 

Remark 1.3. We refer the reader to [20] for theories of stochastic integrals with respect to 
multi-parameter random processes. In fact, the limits above can be written as the following 
Gaussian distributions: 



a ^g<p(x)u (x)dW(x) ul =^M(0,a 2 \\u g^\\ 2 L2 ), (1.15) 
g V ( X )uo(x)W a (dx) distribution f k{v*G V ) g fegg) ^ ). (L16) 

7 y2 \x-y\ a 

Comparing these results with Theorem 12. II below recalling the theory of random fluctu- 
ations in the continuous setting, and with the paths in Fig. [21 we find in the LRC setting 
that the multi-scale scheme (jl.8p captures the theoretical Gaussian limit fluctuations after 
e and h are successively sent to zero. Furthermore, the scheme is robust in the sense that 
it provides the correct fluctuations for arbitrary small patches with < 5 < h (both being 
independent of e). For SRC medium, however, the correct limit for the random fluctuations 
is captured only when 5 = h, that is Kg = K for all K £ 7~h- The amplification effect in 
the case of S < h is again characterized by (h/S)*. The main results hence generalize the 
findings of [8] to a higher dimensional setting. 



7 



Remark 1.4. The main results are stated under the assumptions in Remark 11.11 When 
the ratios |/|X,5|} are not uniform, the limit in ()1.12|) does not have a simple form and 
must account for the non-uniform amplification factors over different triangulation elements. 
Nevertheless, the main conclusions in the above result are not modified. This remark applies 
to the ODE setting in [8] also. 

The rest of this paper is devoted to a proof of the main theorem. Preliminary material 
on random fields and the corrector theory in the continuous scale are provided in Section 
[2j Then main ingredient of the proof is a conservative structure of the stiffness matrix 
associated to the multi-scale scheme; this is considered in section [3j Similar structures 
have been observed and explored in other settings |181 [H] . It allows us to write the discrete 
corrector in the form of oscillatory random integrals. Their limiting distributions are then 
characterized using well established techniques in [TBI 13 E] • This is done in Section HJ 
These sections also include some useful results on the scheme, such as the H 1 estimate of 
the solution to fjl .8[) . which are interesting in their own right. We conclude this introduction 
by several comments. 

1.4 Further Discussions 

This paper studies the specific multi-scale scheme (11 .8f) for the elliptic equation (11. 5ft with 
a random potential. The analysis takes advantage of the conservative structure of the stiff- 
ness matrix. We refer to Proposition 13.41 below for a detailed statement. Other schemes 
possessing this property can be analyzed similarly. To simplify the presentation, we consid- 
ered first-order nodal basis on a uniform triangulation. For higher order schemes in which 
basis functions occupy larger sub-domain of Y, and for general regular triangulation where 
different nodal basis may occupy different number of triangles, the structure in the stiff- 
ness matrix is more complicated. Nevertheless, we believe that the analysis should extend 
without major differences to this more general setting. 

The scheme (|1.8[) fits within the framework of HMM. We refer to [15] for references on 
this method applied to operators of the form C £ = r=i ^a{0'aii{ x /^)dp)-, which are more 
complex to analyze than (jl.5j) . In our scheme, the macro-solver is essentially a standard 
finite element method for the homogenized equation fjl .6[) . and the micro-solver essentially 
estimates the effect of the random potential on small patches Kg- 

Other multiscale schemes and methodologies have been developed for C £ using properties 
of the medium such as separation of scales, periodicity, or ergodicity, e.g. [3j HI [18]. For 
instance, the Multiscale Finite Element Method (MsFEM) in [18] constructs oscillatory 
bases by solving £ e -problems on the supports of the nodal bases {</>^} and uses the so- 
called over-sampling strategy to diminish the resonance errors introduced by the artificial 
boundary conditions of the local £ £ -problems. It would be interesting to investigate how 
random fluctuation are captured by this scheme and in particular what is the effect of the 
over-sampling strategy. The differential operator in (jl.6p does not exhibit such resonances, 
and hence this paper does not address such issues. 

Other multi-scale schemes approach differential operators with rough coefficients like C e 
without assuming any separation of scales or special properties. For instance, [24] constructs 
oscillatory bases by solving £ e -problems on sub-domains that are larger than the supports 
of {(ft 1 - 1 } but still small compared to the whole domain Y . It was proved there, using the so- 



S 



called transfer property of the divergence operator [10], that the resulting finite dimensional 
space can be used to solve the whole £ e -problem with errors that are independent of the 
regularity of {a a p}. Analyzing the fluctuations in such schemes is beyond the scope of this 
paper. 



2 Review of Corrector Theory in the Continuous Scale 

In this section, we review the corrector theories for (|1.5p developed in [161 [5]. They are 
formulated for the following random fields. 



2.1 Random field settings 



In the elliptic equation (jl.5p . the heterogeneous potential, denoted by q £ {x) henceforth, 
consists of a slowly varying part qo(%) and a highly oscillating part q e (x). The latter is 
modeled as g(|,cj), that is, spatially rescaled from some random field q(x,uj) defined on the 
probability space (£l,F,F). In the sequel, E denotes the mathematical expectation with 
respect to the probability measure P. 

We assume that q(x,ca) is stationary. That is to say, for any positive integer k and 
/c-tuple (x\, ■ ■ ■ , Xfc), for any point z and any Borel measurable set icK k , one has 

P{(g(zi), • • • , q{x k )) eA}= + *),•••> V(xk + z)) e A}. 

With this assumption, q admits an (auto-)correlation function R(x) defined by 

R(x) := Eq(y)q(y + x) = Eq(0)q(x). (2.1) 

It is easy to check that R is symmetric and positive definite. Due to Bochner's theorem |27| . 
the Fourier transform of R is a positive Radon measure. In particular, when R is integrable, 
one can define 

a 2 : = [ R(x)dx, (2.2) 
JR d 

and it is a finite non-negative number. Without loss of generality, we also assume that q is 
mean-zero. 

A key parameter of the random field that will determine different limiting correctors is 
the de-correlation rate. It is an indicator of how fast (with respect to distance) the random 
field becomes independent. 

Recall that a random field q(x,cj) is said to be p-mixing with mixing coefficient p if 
there exists some function p(r), which maps 1R+ to K+ and vanishes as r tends to infinity, 
so that for any Borel sets A, B C M d , the sub-cr-algebras Fa and Fb generated by the 
process restricted on A and B respectively de-correlate rapidly as follows: 



sup 



E 0/ - E£ Er] 



(Var f Var tj) 1 / 2 



<p(d(A,B)). 



(2.3) 



Here d(A,B) is the distance between the sets A and B. The function p characterizes the 
decay of the dependence of the random field at different places. We refer the reader to [13] 
for more information on mixing properties of random fields. 
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We consider two settings of random fields. In the first case, we say that q(x, uj) is short 
range correlated (SRC). This means 

(51) q is p-mixing with mixing coefficient p(r) such that p(|x|) G L 1 (K d ). 

(52) \q(x)\ < C so that q e (x) is positive for a.e. uj G Vt. 

(53) In this case, the correlation function R(x) is integrable over M. d and we assume that 
a defined in (|2.2p is positive. 

In the second case, we say that q(x, uj) is long range correlated (LRC). In fact, we consider 
the very specific setting as follows. 

(LI) q{x) has the form $ o g(x), where g{x,ui) is a centered stationary Gaussian random 
field with unit variance and heavy tail, i.e. 

R g (x) := E{g(y)g(y + x)} ~ K g \x\~ a as \x\ -)■ oo, (2.4) 

for some positive constant K g and some real number a < d. 

(L2) The function is uniformly bounded so that q £ {x) is positive for a.e. uj. Further, we 
assume the Fourier transform $ satisfies that f |^>|(1 + |£| 3 ) is finite. 

(L3) We assume also that # has Hermite rank one, that is 

/ <f>(s)e-^ds = 0, Vi := / s^{s)e~^ ds ^ 0. (2.5) 

As a consequence k := V^Kg defines a positive number. For more information on the 
Hermite rank, we refer the reader to |29|. 



2.2 Corrector theory in the continuous scale 

The corrector theory for the elliptic equation with random potential, that is the limiting 
distribution of the difference between u e and uq which solve (11. 5ft and (11 .6f) respectively, 
has been investigated in [161 E] in the SRC setting, and in [6] in the LRC setting. Using 
the notations and random field settings introduced above, the results in dimension two of 
these references can be summarized as follows. 

Theorem 2.1 (|16} [5| [6]). Let u £ and uq be as above and let the dimension d = 2. Denote 
by G(x,y) be the fundamental solution to the Dirichlet problem (jl.6D . When the random 
potential q(x) satisfies the SRC setting, we have 

U e {x) — Uq{x) distribution / „/ \ / \7ttt/ \ /„ m 

> o- G(x,y)u (y)dW(y) (2.6) 



weakly in the spatial variable. When the random potential satisfies the LRC setting, we have 

U e {X) — Uq{X) distribution /— I „, \ / UI , Q / , \ 

'k \ G(x,y)u (y)W a (dy) (2.7) 



\Je a e^o J Y 

weakly in the spatial variable. 

Here, W and W a are the same as in Theorem 11.21 The convergences above are weakly 
in the spatial variable in the sense of (|1.15j) and (11.161) . 
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3 Analysis of the Discrete Equation 



In this section, we analyze the heterogeneous multi-scale scheme (jl.8p in detail. In partic- 
ular, we prove that the scheme with e 5 < h admits a unique solution in the space V h 
that approximates uq in H 1 . With the standard uniform triangulation, we show that the 
stiffness matrix associated to the scheme has some conservative form, which allows us to 
write the discrete corrector conveniently in terms of their coordinates. In the next section, 
we use this discrete representation to prove the main theorem. 



3.1 Well-posedness of the scheme 



The multi-scale scheme (jl.8p with 5 = h coincides with the standard FEM and is well-posed. 
For the sake of completeness, we show that this holds also for 5 < h. 

Recall that V h is the finite dimensional subspace of Hq(Y) with nodal basis {</> y } defined 
in section 11.21 We have defined three quadratic forms: A £ for the heterogeneous equation 
(|1.5p . A £ ,S for the heterogeneous multi-scale scheme which is an approximation of A £ by 



local integration, and A ' which is like A e ' but uses the mean coefficient qo only. In fact, 
Aq' 6 may also be seen as a heterogeneous scheme for the homogenized equation (jl.6p . which 
in turn is associated with the quadratic form 

Aq(u,v) = J Vu-Vv + q uv dx, u,v <E Hq(Y). (3.1) 

Let xk denote the barycenter of the element K £ Th- Exploring the expression of A £ ,S and 
Aq' S , we further define formally 

A%' S (u h , v h )[x K ] = / Vu'" • Vv h + (q + q £ )u h v h dx, 

and similarly define Aq' S (u h ,v h )[xx]- Hereafter, the integral symbol with a dash in the 
middle denotes the averaged integral. 

We define the error due to the heterogeneous choice of integrating element (K$ instead 
of K) by 

e(HMS):= m ax sup ^ $M f ^ ^11. (3.2) 

KeT h V^3uW0 \\ u \\H' I -{K)\\'" \\m(K) 

With this notation we have the following theorem. 

Theorem 3.1. Assume that qo is a nonnegative constant, and that q £ (x) + qo is uniformly 
bounded and nonnegative. There exist unique solutions in V h for the schemes (II. 8p and 
(jl.lOp . Let u £ ,S , Uq ,S be the solutions. Let uq solves (jl.6p . For e <C 5 < h, we have 



\\u h / - uqIIhi < C(h + e(HMS)), (3.3) 
The above estimates hold also if we replace u £ ' S by Uq' S and delete the term e(HMS). 
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Proof. Let p be either e or 0. The existence and uniqueness follow from 

A h /{v h ,v h ) > CII^H^!, for any v h G V h . 
Indeed, because \/v h is constant on K € Th and qo + q £ is non-negative, we have 



A h /{v\v h )> V \K\-f \Vv h \ 2 dx = V / \Vv h \ 2 dx = \v h \ 2 Hl >C\\v h \\ 2 Hl . 



KeT h " Ks KeT h 

Here and in the sequel, | • |#i and | • \ W k, P are the standard semi- norms of the corresponding 
Sobolev spaces. 

We apply the first Strang lemma (Theorem 4.1.1 of P2j), and obtain 

II M|| <r • f /„ fc„ \A h e \v\w h )-A,{v\w h )\ 

\\uq~u' \\ H i<C ml kto — « m + sup n — rj. 

v h eV h V w h eV h \\ w Wh 1 

Set v h = -/Two) the projection of no to the space V h . From classical interpolation result, e.g. 
Theorem 3.1.6 of [12], we have 



||i7«o - uq\\ h i < Ch\\u \\ H 2. 

For any w h 6 V' 1 , we have 

\A h /{v h ,w h ) - A (v h ,w h )\ < \A h /{v h ,w h )-A h /{v h ,w h )\ + \A^ s {v h ,w h ) - A {v h ,w h )\. 

For the first term, we have 

\A h /(v h ,w h )-A h /(v h ,w h )\ < Y, \K\\A h /{v\w h )[x K ]-A h \v\w h )[x K ]\ 

KeT h 

<e(HMS) Y WAhHkMWhHk) 
KeT h 

< e(HMS)||u fc || ff i||«; h ||Hi- 

In the equalities above, we used the definition of e(HMS) and Cauchy-Schwarz respectively. 
For the second term, we first observe that 



A h /(v\w h )-Ao(v\w h )= Y\w\{I 



q v h w h dx-\K s \(q v h w h )(x K 

I K S 



q v h w h dx - \K\(q v h w h )(x K ) 

K 



The items in the sum can be recognized as errors of barycenter numerical approximation of 
integrals. Error estimate for such numerical quadrature is discussed in the next lemma and 
by (pT3j) we have that \Aq ,5 {v h , w h ) - A {v h ,w h )\ is bounded by 



V) 

,2 



^K^l^W^iK) 



Yl C ho\\d \ pS\\v h \\ H i (Ks) \\w h \\ L 2 {Ks) + h\\v h \\ H 
KeT h 1 

<C'h\\q \\ c i Y \\v h \\m{K)\\w h \\ L 2(K) < Ch\\qo\\c i \\v h \\m\\w h \\m- 

Ken 
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Combining the above estimates, we find that 

IN - Ue' 5 \\m < \\nu - u \\m + (e(HMS) + Ch)\\n Uo \\ H i < C(h + e(HMS)). 

The constant depends on H^ollc 1 ) 1 1 u o I Iff 1 an d some uniform bound of h/5 and hence is 
independent of e, h or 5. □ 

The following lemma concerns error estimate for barycenter numerical quadrature of 
product of two functions in Pi {K) , the space of linear polynomials on a triangular element 
K. It is stated in the simplest setting thought it can be generalized to regular element 
easily. This lemma is used in the proof of the previous theorem. 

Lemma 3.2. Let K be an isosceles right triangle with unit side length. Let K be the image 
of K under the linear transform F{x) = Bx + b G M 2 . For any function tp on K, let ip be 
the function ip o F. Assume qo G W l,0O (K). Then for any v,w G Pi(K), we have 



qo(x)v(x)w(x) dx — \K\ (qovw)(xK) 



< cil-SlllkollvK 1 ' 00 ^)!^!!//^^)!!^!^ 2 ^)- (3-4) 



IK 

Here, xk is the barycenter of K; \\B\\ is the matrix norm of B. 

Proof. We follow the steps in the proof of [121 Theorem 4.1.4 ]. Consider any ijj G W 1,oc (K) 
so that ip = ipoF is in W^°°{K). Let \E(il)w)\ denote the error of the barycenter quadrature 
for the integral J K ipwdx. After change of variables, 

E(ifnv) = |det(J3)| (J 4>(x)w(x) dx - \K\(^w)(x^)\ = \det(B)\E($w). 

On the reference element K, since all norms on Pi(K) are equivalent, we have 

\E(i>w)\ < C0\\ Loo{&) \\w\\ LOO{&) < C\\i>\\ wl ^ {k) \\w\\ L2{ ^ y 

We view E(- w) : tp i— )• E(t/jw) as a linear functional on W 1,QO (K). The above estimate 
shows that E{- w) is continuous with norm less than CH?/)!^^)- We check also that E{- w) 

vanishes on Pq(K), the space of constant functions on K. Therefore, due to Bramble-Hilbert 
lemma \12\ Theorem 4.1.3], there exists some C such that for all ijj G W 1,00 (K), 

\E(^w)\ < C\\E(- *)|| £(H /i,cx 1( i- )) |V'| H/ i, o ( ^ ) < C\\w\\ L2{k) \Tp\ wUoc{k) . 
Take ip = qov. We check that 

\lp\ W l,oo^) < l?olw 1 ,°°(K)Flll,°°(K) + 11^0 11^00(^)1^1^1,00^ 
- C (l90| H /l,oo(^-)||«|| L 2(X) + ll90|| L oo(^)|w| H 

The last inequality holds because v G P\{K) and all norms on P\{K) are equivalent. Finally, 
recall the relations [12} Theorem 3.1.2] that for any integer m > 0, any q G [l,oo], and for 
any ip G W m ' p {K), 

_i 

\ip\ Wm , g{ K) < C\\B\\ m \det(B)\ *\v\ W m, q{K) . 
Combine all the estimates above, we obtain the desired estimate. □ 
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For the heterogeneous multi-scale error, we have the following result. We do not intend 
to make these estimates sharp. 

Theorem 3.3. Let dimension d = 2 and e <C 5 < h. Let e(HMS) be the multi-scale 
heterogeneous error defined in (|3,2|) . we have, for e > sufficiently small, the following 
estimate: f 

C ( - j 5 , in the SRC setting, 
E e(HMS) < <J yd \ (3.5) 

C (^j 2 , in the LRC setting. 
The constants C above does not depend on h,5 or e. 

Proof. For any v h 6 V , its restriction v h \k is in P\(K). Therefore, according to the 
definition of e(HMS) in (j3.2|) . we try to estimate 

\K\\A h i\u h ,v h )[x K ]-A h Q \u\v h )[x K ]\ 

eK '■= SUP 7 7 7 7 . 

v,weP\{K) \\ v \\h 1 (k)\\ w \\h 1 (k) 
Since Ae' S — Aq' S may be viewed as a bilinear form, we have 

3 1 

e K <[ \K\MV^rn^n)[x K ]-A h '\(t> m An)[xK]t 
m,n=l 

where {(fi m ,m = 1,2,3} are some orthonormal basis of P\{K). Let X^ nn denotes the 
difference \K\ (Ae' 5 ((l) m , <j)n)[%K\ ~ A^^m, M[xk]) ■ By Cauchy-Schwarz, 



3 

£ |2 

m,n 

m,n=l 



Eejf<[ ^2 E\X. 



Therefore, we estimate E|X^ n | . From the definition, we check that 



X mn= q e (x)lpm,n(x)dx, with ^ m ,n{x) = %K S {x)}—r(j) m {x)(j) n {x) . 



\K\ 



Here and in the sequel, \A is the indication function for a set icK 2 . Abusing notations, 
we use X £ and ip instead of X^ n and ^ m ,n in the following. 

Let us estimate E|X e | 2 . In the SRC setting, it is known that e~^X e converges in 
distribution to a mean-zero Gaussian variable with variance a 2 see [5j Theorem 3.8]. 
Therefore, for sufficiently small e, we have 

E\X e \ 2 < Ce d \\R\\ L 4l>\\ 2 L2 = C^\\R\\ L ^) M \\<f> m Mh{K t ) < C\\R\\ L ,h 2d ( £ -) d . (3.6) 

Here R is the correlation function of q defined in (|2.ip . We used the fact that l-R^/l-K^I = 
(h/S) d and we calculated that ll^Vre^nll 2 ^^) 

< \K 5 \ < C8 d . 
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In the LRC setting, it is known that e~ 2 X £ converges in distribution to a mean-zero 
Gaussian variable with variance <g> ip\\i J 1 (Y 2 ,K,\x-y\- c 'dxdy)'i see El Lemma 4.3]. Conse- 
quently, for sufficiently small e, we have 



K 5 xKs \x-y\ 

Ce a (j) 2d U m cl )n \\ 2 4 <Ch 2d (~) a . 

L 53 " 



(3.7) 



In the second inequality we used Hardy-Littlewood-Sobolev inequality \22\ Theorem 4.3], 

A — a. 

and we calculated that I m n ] 4 < liT^I^^. 

The inequalities fl3l| and (^7) show that E e K is of order /i d (f) d / 2 and h d (^) a / 2 in the 
SRC and LRC settings respectively. From the definition of e(HMS), we see that 

E e(HMS) < V E e K < sup E e K . (3.8) 

Here, p- is a bound for the total number of elements in Th- Since the estimates (|3.6p and 
(|3.7p are uniform over K £ Th, we obtain the desired estimates. □ 



3.2 Coordinate representation and conservative form 

The next step is to reformulate the multi-scale schemes (|1.8|) and (|1.1U|) as linear systems for 
the coordinates of the solutions in V h , to investigate the structure of the associated stiffness 
matrices, and to write the discrete corrector Us' S — Uq' S in terms of their coordinates. 

We start by introducing some useful notation. In the triangulation illustrated by Fig. [3l 
we identify each grid point (ih,jh) with a unique two dimensional index The set of 

inner grid points are denoted by I = | 1 < i,j < N — 1}, and the set of all grid 

points including the boundary ones is denoted by I = | < i,j < N}. We define six 

difference operators df : I — > I as follows: 

df(i,j) = (i±l,j), d±(i,j) = (i,j±l), df(i,j) = (i±l,j±l). (3.9) 

Here, s = 1,2,3 denotes three directions: horizontal, vertical and diagonal; the plus or 
minus sign indicates forward or backward differences. 

In the sequel, we often write (i, j) simply as ij. For each ij £ I, there corresponds a basis 
function (j) 13 which is piecewise linear on each element K £ Th, has value one at ij and has 
value zero at other nodal points. Any function v h in the space V h can be uniquely written 
as v h {x) = Vij < i >l ^( x )i an d the vector (Vij) G ]R( Ar ~ 1 ) x ( 7V ~ 1 ) is called the coordinates 

of v h . We identify R^- 1 )*^" 1 ), the space for the coordinates, with V h itself. Now, the 
difference operators induce difference operators Djr on V h as follows: 

ikv.j v dfij v, r djv* v u v LU . (3.10) 

Note when d^ij lands outside of X. i.e. on the boundary, the value V A ±-- is set to zero. 
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Using the coordinate representation of functions Us ,S = £V . UfjCp^ and Uq' S = U^cpi , 
we can recast the heterogeneous multi-scale schemes (jl.8p and (ll.lOp as the following sys- 
tems: 

Af jkl U £ kl = (f,^), (3.11) 
A° ijkl U° kl = (f,^). (3.12) 

Here, the stiffness matrices are defined by 

A! jkl = A h /{^^% A% kl = A h /{^A kl )- 
These stiffness matrices have the following structures. 

Proposition 3.4. Let A p = (A p - kl ) with p = or e be the stiffness matrices above. We 
observe 



{ij}U{dfij | s = 1,2,3}. 

8=1 

for some of- that can be explicitly computed as in f|3. 14j) below. 

Proof. The first two observations are obvious, so only the third one needs to be stressed. Ac- 
cording to (jl.8p and (jl.lOp . to calculate A p ^- we need to integrate the function |V^ j (:e)| 2 + 
q p {x)\4>^ {x)\ 2 . We observe that the support of (/>* J , denoted by /Qj, is a hexagon consist- 
ing of six triangle elements as illustrated in Fig. 0}-Left. The integration is actually taken 



(PI) A ijkl ~ A kiij>' 
(P2) A p m = unless kl G := 
(P3) For any ij S I, we have 
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over K s ip the region obtained by shrinking the triangle elements in /Qj with respect to 
their barycenters as illustrated in Fig. 0]- Right. Let us consider a typical triangle in K-ij 
with nodal points ij, d^ij and dtij. Abusing notation, we call it K and the corresponding 
smaller triangle Kg. Note Kg corresponds to the shaded region in the figure. On this region, 
the three non-zero basis functions are cf) lJ , </> d i ij and 4> d 3 ij . They satisfy 

+ cp d i ij + (p d z ij = 1, V<p ij + Vcp d i ij + V4> d t^ = 0. 

Multiply <f> %3 on both sides of the first equation, and V on the second equation. We have 

(0«) 2 = pi - (tj) d tv + iftfV)^, \Vf j \ 2 = -(V0 d ** j + V) ■ 

Consequently, we have 

A%' 6 ((t> ij ,<t> ij )[x K ] = / |V^| 2 + q p \^\ 2 dx 



• p (j) ij dx - V / Vcf) ij • Vcf) dtij + q p (j) ij <j) dtij dx 



s=l,3" 

<fct> ij dx - A h /{^^ Atij )[x K ] - Ai' 6 (f j A dtij )[x K ]. 



K 6 

Summing over the integrals on all six triangles, and using the notations of A p , Ap' S and 
Ap' 6 , p = 0, s, we see that (13.131) holds with d v - defined by 



V \K\-f q p <j) ij dx. (3.14) 



4 

'K S 



This completes the proof. □ 
It follows immediately that the matrix AP acts on vectors in V h as follows: 

3 



(A p V)ij =Y J Dt{*lfD-V ij ) + d%V ih 



8=1 



where a- f is short-hand notation for A p _ and it has the expression 

l 3 ijd s ij 



:= J2 \ K \-f Vo" • Vo' L '' • </'VV ; * ij dx. 



s,p 

<*ij 

'K S 



Note that when dfij lands outside of I, i.e. on the boundary, (j) ds %3 is the unique continuous 
function which is linear on each K G Th, has value one at d^ij and value zero at all other 
nodal points. Finally, taking the difference of A 6 and ^4° we obtain 

3 

{A e V _ A o v) .. = £ D+^D-Vij) + d £ij V ij: (3.15) 

8=1 
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where the vectors (ct £i j) and (d e ij) are 

a s Eij :=atf-a°f= £ K /' q e ^^dx, (3.16) 

d £ij := dfj - d° tj = £ \ K \-f Qefldx. (3.17) 

Formula f)3. 15j) is essential in our analysis because it provides an explicit expression of 
the discrete corrector u e ' — u Q ! . Identify these solutions with the vectors (JJfj) and {Ufj) 
in (|3.1HI3.12p . We verify that 

A% Jkl (U e -U°)ki = -(A e -A ) ijkl Ut l . 

Here and after, we will use the summation convention that repeated indices are summed 
over unless otherwise stated. Let (G?- fc j) be the inverse of ^4°. We then have 

(C/ £ - U°)ij = -G% kl {M - A°) klst Ul t . (3.18) 
Using the formula (|3.15p and summation by parts, we obtain 

(U* - U% = -G m {A s - A°) klts U? s - G ljkl (A £ - A°) kUs {W - U°) ts 

3 

= J2( D sG ijk i)(a s £ D-U°) kl - G m (d £ U°) kl (3.19) 

s=l 
3 



+ ^2(DjG ijkl )(a s £ Dj(U £ - U°)) kl - G ljkl (d £ (U £ - U )) 



ki- 



Note that for two vectors of the same dimension, say d e and U° above, the notation (d £ U°) 
is the vector obtained by multiplying the corresponding components. This decomposition 
formula will be the starting point of our analysis in the next section. 



4 Proof of the Main Results 

In this section, we prove Theorem 11.21 using the coordinate representation (|3.19p of the 
discrete corrector. 

From (|3.16p and ()3.17j) we see that a £ and d £ may be seen as averages of fluctuations 
and hence are asymptotically small. In the decomposition formula (|3.19|) . the first sum 
involves linear terms of these random processes, while the second sum involves product of 
a £ or d £ with U e — U°. The second term hence is much smaller if we can control U £ - U°. 
This is done in the following lemma. 

Lemma 4.1. Let £/?„• and £/£ be the coordinates of the numerical solutions to the random 
% j 1 j J 

and the deterministic equations (I1.5P and (I1.6P respectively. Suppose that there exist some 
constants C > 0,jj E = 1, ■ ■ ■ ,4 so that 

\D7G m \ < ChPADjUfA < Ch^,\G ijkl \ < Ch^and {UfA < Ch 14 (4.1) 
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for any s = 1, 2, 3. Let d = 2. T/ie following holds. 

(1) // i/te random process q satisfies the SRC setting, we have 

E\\U £ - U°\\% < Ch 2( - min{ll+ ^ 2 ' l3+li}) \\R\\i (4.2) 

(2) If the random process q satisfies the LRC setting, we have 

E\\U £ - U°\\] 2 < C(a, K )/ l 2 ( min {7i+72,73 +74}) f^y _ (4.3) 

The constant C does not depend on h,5 or e. 
Proof. We observe that (|3.19p can also be written as 



3 



(U £ - U%j = Y,(D;G i:i ) kl (a s £ D;U £ ) kl - G ljkl (d £ U £ ) kl . 

s=l 

Using the bounds in (|4.1|) and Cauchy-Schwarz, we have 

3 

E\U £ -U ^ < CN 2 h 2 ^ 1+ ^^^E\a s £kl \ 2 + CN 2 h 2 ^ + ^^2E\d £k i\ 2 . (4.4) 

s=l M kl 

Here, iV 2 is the number of nodal points, i.e. N 2 = \X\ ~ h~ 2 . Take expectation. It suffices 
to estimate E\a s £kl \ 2 and E\d ek i\ 2 . We rewrite (|3.16«3.17p as 

®M = J Qe(x)a s kl (x)dx, d ekt = J q e (x)b k i(x)dx, (4.5) 

with a s e and b e defined by 

4i(*)= E l§- { XK s (x)^ l (x)^ kl (x), b kl (x)= £ -H x Kg (x)4J*(x). (4.6) 

Note that |i^|/|ET,5| = (h/5) d . Note also that a k , and b k i are uniformly bounded on Y. 
Hence, we recognize a s £kl and d ek i as oscillatory integrals of uniformly bounded functions 
against fast varying mean-zero random processes. Such integrals are well understood. In 
fact, a s E has the same form as X £ in the proof of Theorem 13.31 and can be estimated in the 
same manner. In the SRC setting, we have that 

E|c4,| 2 < Ce d \\R\\ L 4a s kl f L2 < C\\R\\ L ih 2d ( £ -) d . (4.7) 

In the LRC setting, the above estimate should be changed to 

E\a s £kl \ 2 < Ce a \\a s kl ® a s kl \\ L i {YxY>lx _ yl ^ dxdy) < C(a, K)h 2d (^) a . (4.8) 



The mean square of d ek i can be similarly estimated. Substitute these estimates into (14. 4j) 
to control the mean square of (U £ - U°)ij; note that the sum over kl introduces a factor 
of hT d which is the number of items in the sum. The estimates of (U £ - U%j are uniform 
in ij, summation over ij yields the desired results. Note that this additional summation 
introduces another h~ d to the estimates. □ 
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Lemma 4.2. Under the same condition of the previous lemma, we have 

3 

(U e - U% = J2(DjG ijkl )(atDjU°) kl + G ljkl (d £ U°) kl + 4. (4.9) 

8=1 

Further, the error term rfj satisfies 



C/l mm{ 7l , 73 }+min{ 7l + 72 , 73 +74 }+§ ^ d ^ ^ ^ 

C/l mm{ 7l , 73 }+mm{ 71 + 72 , 73 + 7 4}+| Q° j m ^ g £J j C settin5 . 



Proof. The decomposition holds with 

3 

4 = E E(^~ G ^a^£>8~(£4 £ * " ^% + E G mdeki(U £ - U°) kl . (4.11) 

8=1 kl kl 

Bound the D~Gij k i and Gij k i terms by (|4.1|) . and use Cauchy-Schwarz. We get 

3 

141 < ch? 1 (u £ - u°)\\p + ch ri \\de\\ P \\u £ - u°\y. 

8=1 

Note that \\D~ (U £ - U°)\\j 2 < C\\U £ - U°\\ £ 2. Take expectation and use Cauchy-Schwarz 
again to get 

3 

E|4| < Ch^ m\ a eWh E \\U £ - U°\\%Y +Ch^ (E||d £ |ll2E||C/ £ - U°\\ 2 e2 ) 2 • (4.12) 
8=1 

Summing over kl in the estimates (|4.7|) and (|4.8|) . we have 

f C\\R\\ L ih d ^-) d , in the SRC setting, 

n<\\l < I e 

C{a, K)h d (-) a , in the LRC setting. 
v 

The same estimates hold also for E||d e ||| 2 - Substituting these estimates, together with 
(I4.2M.3I) . into (|4TT2l) completes the proof. □ 

Remark 4.3. The assumption (|4.ip is not a restriction because jj there can be negative. 
Indeed, consider an arbitrary triangle element K. Without loss of generality, let its vertices 
be {ij, i — + 1}. For any function v h E V h with coordinate vector V, we verify that 

Cih 2 (\Vij\ 2 + l^-i/ + \V ij+1 \ 2 ) < \\v h \\ 2 LHK) < C 2 h 2 (\V lJ \ 2 + IV-^I 2 + \V lJ+1 \ 2 ). (4.13) 

This is due to the fact that all norms on the finite dimensional space V \k are equiv- 
alent and h? is the right scaling. We can check also that \7v h \K is a constant vector 
(D^Vij,D^Vij^i)/h. It follows that 



\V\w^i{K) 



Ch<- 1 \\{DiV ij ,DtVi j+ i)\\ q = CfcHflDJ-W + \D2V ij+ i\ q y. (4.14) 
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Now for lie' 6 , we know its H 1 norm is bounded independent of h and e. Applying the 

results above we find that \UfA < Ch~ l and \D7~UfA < C. Other elements of U e and D~U £ 

j 1 1 L j 1 a 

can be estimated in the same way. Hence, we may choose 72 = and 74 = —1. Similarly, 
the discrete Green's function G h {x,y) = Gijki^ \x)(f> kl (y) is known to have W 1,q norm for 
some q < 2 bounded by C|log/t| for any fixed x; see [171 Theorem 5.1]. Using (|4.13j) and 
(|4.14p we may choose 71 and 73 properly, say 71 = 73 = — 1 • □ 

Now we are ready to prove the main theorem. We denote the solution operator of the 
averae ed scheme (fLTU]) by G h ' 8 . That is, G h > 5 f is the unique solution in V h that satisfies 
(jl.lOp . Using the coordinate representation and summation convention, we verify that 



g h ' 5 f(x) = G ijM {f,(j) kl )<i) v {x) 



kl\ 



(4.15) 



h,6 



Proof of Theorem \1.2[ Take any test function (p £ C 2 (Y). Let us denote the function Q 
by m h and its coordinate vector by (M^). Let /3 = d in the SRC setting and (3 = a in the 
LRC setting. Prom the decomposition (14. 9p . we write 



1 



UP 
1 

w 

1 



y{x)[u h /-u n Q '°}dx 



1 



:(U £ 



^2D7G ljkl (a s £ D7U°) kl + G ijkl (d E U°) kl + 



\s=l 



V DjM kl {a s £ DjU°) kl + M kl {d £ U G 



s=l 



ikl 




(4.16) 



In the last equality, we used the fact that Gij k \ = G k uj and recognized the coordinate M k \ 
according to the formula (I4.15p . 

First convergence as e — > while h is fixed. Let us control the last term first. Thanks 
to the estimate (|4.10p . we have 



E 



1 



TP 11 



< 



iE\r!A)\( V ,^)\<C(h)y\\ L We^ 



(4.17) 



This term hence converges to zero in L X (P) and does not contribute to the limiting dis- 
tribution. So we focus on the other two terms that are linear in a s £ and d £ , respectively. 
Substituting the oscillatory integral representations (|4.6p into (|4.16p . we find that 



y 



h,S 

£ 



h,S 



]dx 



7= / q £ {x)L^ S (x)dx + -= / q £ {x)L 
Up Jy vep Jy 

q £ (x)L h ' S (x)dx. 



2' S (x)dx 



(4.18) 
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Here, L h - ,& , j 



1,2 and L 



h.5 



L^' S + Lg' 5 depend on cp through M and are defined by 



^2D;M kl (a s kl (x)D;U c 



>ki 



s=l 



E E ]§iXK s (x)J2(DJM kl )(DJU^ 

kl K£lC kl 1 51 8=1 



L 2 S ( X ) 



bki( X )(Mu») kl =Y: e tSi 

fcZ K&C kl 1 d| 



(4.19) 



X ^(x)M H C/ fc V fe '(x) 
K 



v IS 

j4 1*1 



^ I- 



X ^(x)77' l (mV h ' 5 )(x) 



Here, Xk contains the indices so that (kh,lh) is in K and TI h {m h u^' 5 ) is the projection 
in V h of the function m h UQ' 5 . Now the convergence in item one and the first conclusion 
of item three of Theorem (|1.2p follows from the representation (|4.18p and the well-known 
results on limiting distribution of oscillatory integrals; we refer the reader to Theorem 3.8 
of [5] for the SRC setting, and to Lemma 4.3 of [6] for the LRC setting. 

Second convergence as h — >■ 0, SRC setting. Now we prove item two of the theorem. 
It concerns the limiting distribution, as h goes to zero, of the Gaussian random variable 
which is obtained as the limiting distribution in the first step. This step depends on the 
correlation length of the random field and needs to be considered separately for the SRC 
and LRC setting. Here we focus on proving (11.120 first. 

We have the following key observation: 

LY — > in L°°(Y) as h -> 0. (4.20) 
Indeed, for any fixed x G Y, since \(p l ^\ < 1 uniformly and l-R^/l-fQl = (M -1 ) 2 , we have 



\L h /(x)\ < C 



h 



3 

2 E 

s=l 



D~M kl 



h 



DsU° kl 



h 



<C - h 2 \m h \ H1 \u h /\ m . 



Since Uq' S and m h are yielded form the scheme (jl.lOp for smooth right hand side / and tp, 
they have bounded H 1 norms. We assume that the ratio h/5 is fixed while h is sent to zero. 



Therefore, the above estimate shows that L^' S goes to zero uniformly, proving the claim 



According to (|4.18|) . the left hand side of (|1.12p can be written as 

a [ L h /{x)dW{x) + a [ L h /{x)dW{x). 
Jy Jy 



(4.21) 



Our plan is to show that the second term above converges to the right hand side of (|1.12p 
while the first term above converges in probability to zero; this is indeed sufficient for (|1.12p . 
Since all random variables involved are Gaussian, we only need to calculate their variances. 
Thanks to Ito's isometry, we have 



Vara J L h l '\x)dW{x) = a 2 j \L h / {x)\ 2 dx . 
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Due to (|4.20p . the above variance goes to zero, proving our claim for the first term. For the 
second one, we have again 



Vara / L^ S (x)dW(x) = a 2 / \L h /{x)\ z dx 



T) 2 £ \K\j K \n\m h uY){x)\ 2 dx. 



Y JY \ o J ^ J Ks 

We recognize the sum in the last term as a barycenter approximation of the integral that 
gives the L 2 norm square of II h (m h UQ ,S ). Thanks to Lemma 14.41 below. \\II h (m h UQ' S )\\i J 2 
converges to ||ito£/y||i,2. This implies that the variance of the second term in (|4.2ip converges 
to (ah/5) 2 \\u Q^\\ 2 L2 , proving (|1.12j) . 

Second convergence as h — >■ 0, LRC setting. Now we prove (j!.14|) . Like in (|4.2ip . we 
can write the left hand side of ()1.14p as a sum of two Gaussian random variables. Using a 
modified isometry, we write the variance of the first variable as 

Var a [ L h /{x)W a {dx) =11 kL \ {x)L } & {v) dxdy = S(L h /). 
Jy ' J Jyz \x-y\ a 

4 

Here, we define the operator : L 4 -<* — y R as 

S(g) := \\g ® g\\Li(Y2, K \x-y\-<*dxdy) = J J \ x _ y \a dx ^' ( 422 ) 

Recalling the Hardy-Littlewood-Sobolev inequality, Theorem 4.3 of [22], we have 

\,y(g)\ <KC{a)\\g\\ 2 , . (4.23) 

Due to (I4.20p . this term goes to zero again and does not contribute to the limiting distri- 
bution. For the contribution of L^' 5 , we have 

h,5 1 \ T h,S I 



Vara| L h /(x)W a (dx) = jj ^ — dxdy 



= T y \K\ 2 -[ J Kn\m h u h > 5 )(x)n h (m h u h > 5 )(y) ^ 

We recognize the last sum as the barycenter approximation of J?(II h (m h UQ ,S )). Now (|4.23p 
shows that J? is continuous on L 3 ^. Since a < 2 and < 2, we have the inclusion 

4 

L 2 (Y) C L A ~ a (Y). Therefore J? is also continuous on L 2 (Y). Applying (I4.24D with fi = ip 
and /2 = /, we conclude that J?(n h (m h UQ ,S )) converges to ^(uoQip). This proves (|1.14p 
and completes the proof of the theorem. □ 

It remains to prove the following key lemma concerning the convergence of product of 
solutions yielded from the averaged heterogeneous multi-scale scheme (jl.lOp . 
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Lemma 4.4. Let Q h ' be the Green's operator of the scheme (jl.lOp . For any two functions 
fj E C 2 (Y), j = 1,2, Zet n h (G h ' S f x Q h ^ / 2 ) 6e toe projection in V h of the product of G h ' S h 
and G h,S f2- We have that 

n h(gh,5 fi gh,6 h) g h g f2t ash ^ . ( 4 .24) 

^4s before, Q above is the Green's operator of the homogenized equation (II. 6p , 

Proof. To simplify notation, let us denote the function Q h,s fj by the functions Qfj by 
itj, j = 1,2. 

The key to the proof relies on L°° error estimates for finite element methods. Such 
results are classic for the scheme with h = 5 as proved in |26t [28] . For 8 < h, as explained 
before we may view the scheme as the standard finite element with (barycenter) numerical 
integrations. L°° error estimates for such practical schemes are more involved but were 
obtained in |30|. 117] , In particular, the piecewise linear FEM with numerical quadrature 
was considered in Theorem 5.1 of [T7j, and it was shown that 

\\uj - Uj\\L°° < Ch 2 \ log/i|||/j||vi/2,oo. 

Since Uj, j = 1, 2, are bounded, the above also implies that 

\\uiit2 - uxu 2 \\l°° < Ch 2 \ log h\\\fj ||^2,oo. (4.25) 
In fact, Theorem 5.1 of |17| also shows that 

Htf/l.oo < WUjWiyl,^ + Ch\l0gh\(\\Uj\\ W 2,ao + \\fj\\ W 2, X ). 

Here, Uj is the FEM solution with h = 5. The above estimate shows that Uj is in W^. 
Since Uj are bounded, we check that u\u\ G W^, . From classical interpolation estimates, 
e.g. taking k = m = 0, p = 00 and q = 2 in Theorem 3.1.6 of [12] . we have 

\\u\u\ - nk(u1u%)\\ LHK ) < C\K\^h\u\ul\ w ^ {Y) . 
Here, IP\ is the projection on the triangle element K. Summing over K E 7~h, we have 

||fiffi§ - n h (u^)\\ L 2 {Y) < Ch\\u1u%\\ w i,oo. (4.26) 
Note that (|4,25p controls Hu^u^H^i^. Sending h to zero, we finish the proof. □ 
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